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Abstract 

In the vicinity of their glass transition, dense colloidal suspensions acquire 
elastic properties over experimental timescales. We investigate the possibil¬ 
ity of a visco-elastic flow instability in curved geometry for such materials. 
To this end, we first present a general strategy extending a first-principles ap¬ 
proach based on projections onto slow variables (so far restricted to strictly 
homogeneous flow) in order to handle inhomogeneities. In particular, we 
separate the advection of the microstructure by the flow, at the origin of 
a fluctuation advection term, from the intrinsic dynamics. On account of 
the complexity of the involved equations, we then opt for a drastic simpli¬ 
fication of the theory, in order to establish its potential to describe insta¬ 
bilities. These very strong approximations lead to a constitutive equation 
of the White-Metzner class, whose parameters are fitted with experimental 
measurements of the macroscopic rheology of a glass-forming colloidal dis¬ 
persion. The model properly accounts for the shear-thinning properties of 
the dispersions, but, owing to the approximations, the description is not fully 
quantitative. Finally, we perform a linear stability analysis of the flow in the 
experimentally relevant cylindrical (Taylor-Couette) geometry and provide 
evidence that shear-thinning strongly stabilises the flow, which can explain 
why visco-elastic instabilities are not observed in dense colloidal suspensions. 
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1. Introduction 


1.1. Observations 

Take a small amount of carbon black powder and disperse it into water: 
this gives pigmented ink. From a rheological perspective, it is a colloidal 
suspension which flows similarly to water, albeit with a somewhat higher 
viscosity. But this Newtonian behaviour, which holds generically for very 
dilute suspensions, is strongly altered when the volume fraction 4> of colloids 
gets larger. Most strikingly, the viscosity and relaxation time of the suspen¬ 
sion increase dramatically when <f> approaches a “critical” packing fraction 
cf>g [(f) g ~ 0.56 for hard-sphere-like colloids [f|). For (j) > <p g , the material 
retains elastic properties over any experimental timescale, in a fashion rem¬ 
iniscent of the emergence of glassiness in supercooled melts of some metallic 
alloys, when the temperature declines [2]. Despite these dramatic changes, 
the structure of the material remains essentially liquid-like throughout the 
transition. Accordingly, it appears sensible to compare the rheology of very 
dense colloidal suspensions to that of other visco-elastic liquids. In particu¬ 
lar, one may wonder why a variety of complex fluids among the latter, such 
as worm-like micelles or polymer solutions [l, 3], are prone to a (non-inertial) 
flow instability in curved geometry, leading for instance to the formation of 
vortices, while, to the best of our knowledge, no such visco-elastic instability 
has ever been reported in very dense colloidal suspensions. 

1.2. A microscopic approach using mode-coupling theory 

The level of difficulty required to rationalise the rheology of suspensions 
strongly depends on the volume fraction 4> of interest. In the dilute regime, 
the fluid is Newtonian; its viscosity is independent of the applied shear rate. 
More quantitatively, the linear corrections to the solvent viscosity due to 
the colloids were worked out by Einstein a little more than a century ago, 
under the assumption of non-interacting colloids |5]. By a detailed study 
of the probability distribution function of particle pairs, the approach was 
extended to interacting colloids by Batchelor and others [6|, and led to a 
description of the semi-dilute regime [7]. For (j) ~ cf> g , collective effects be¬ 
come paramount, in that glassiness can be thought of as the entrapment of 
particles in the “cages” formed by their neighbours; these effects turn a first- 
principles derivation of the macroscopic rheology into a formidable challenge, 
all the more so as the presence of flow distorts the structure of the material 
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away from its “quiescent configuration” and gives rise to complex interplays 


Nevertheless, at the expense of some uncontrolled approximations, the 
mode-couplin g th eory developed by Sjogren, Bengtzelius, Gotze, Sjolander, 


and others [9|, [10] succeeded in rationalising the phenomenology of the glass 
transition by focusing on the evolution and the relaxation of the (slow) den¬ 
sity modes of the system and on their coupling to the other (faster) variables. 


In the last decade, Fuchs, Cates et al. [Ill, |12|. and Miyazaki et al. [13] in a 


parallel endeavour, were able to extend this framework to situations of flow, 
in which the colloids are dragged by a prescribed solvent flow. The state 
of the art of this theory encompasses arbitrary, potentially time-dependent 
incompressible solvent flows, in two or three dimensions |l2l. Il4|. 

However, the derivation hinges on the assumption of a perfectly homo¬ 
geneous flow throughout space; this hampers the investigation of any flow 
instability. Indeed, perturbations, which break homogeneity, are not handled 
adequately; in particular, the mechanism describing their (expected) advec- 
tion with the flow is still missing in the equations. Moreover, the complexity 
of the final equations giving the stress as a function of the strain history is 
a deterrent to any stability analysis in non-trivial geometry. 

1.3. Objectives of the article 

In this contribution, we first propose a general way to extend the formal¬ 
ism and handle flow inhomogeneities, insisting in particular on the recovery 
of a fluctuation advection term and on the limit of locally homogeneous flow. 
Then, we follow the endeavour pioneered in Ref. [15] to reduce the final equa¬ 
tions to a tractable constitutive equation. This will come at the expense of 
very strong (but explicitly exposed) approximations and clearly undermine 
the accuracy of the description. Nevertheless, the ensuing simple model, 
which falls in the White-Metzner class (l6| . will allow us to capture the ex¬ 
perimentally measured low-shear-rate rheology and high-shear-rate rheology 
in a model colloidal glass-forming dispersion 17, 0. Finally, a linear sta¬ 


bility analysis of the flow will be performed, in cylindrical (Taylor-Couette) 
geometry and the (stabilising) effect of effect of shear-thinning on the visco¬ 
elastic flow will be numerically assessed. 


2. Theoretical probabilistic framework 

We start by presenting the theoretical underpinning of the rheological 
equations that extend quiescent mode-coupling theory. 
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Let us consider an assembly of N colloidal particles dispersed in a solvent 
and evolving by Brownian motion in a volume V, for instance with periodic 
boundary conditions. 


2.1. From the overdamped Langevin equation to the Smoluchowski equation 

To describe the microscopic motion of particle i G {1,... ,N}, we posit 
an overdamped Langevin equation acting on its velocity rp. 



solv 



Fi + fi th 


(1) 


Here, Fi is the conservative force that derives from the global potential en¬ 
ergy of the system, and the /j th ’s are random Gaussian thermal fluctuations, 
viz., (fi th (t)) = 0 and (fi th (t) ® fj th (V)) = 2ksT( i 5ij5 (t — t') I, where I is 
the identity matrix in d dimensions. The frictional force on the left-hand side 
(lhs) involves a frictional coefficient ( and the particle velocity relative to the 
(prescribed) local solvent velocity u solv (r); for an incompressible flow, the 
latter should satisfy V • v solv = 0. Hydrodynamic effects, being presumably 
subordinate to short-range interactions in dense systems, are neglected. 

Rather than focusing on the motion of individual particles, we adopt a 
statistical approach. Equation Q] is recast into the following equations for 
the evolution of the probability -0 (T; t) to find the system in the microscopic 
configuration T = {r\,.. . , rjy) at time t [IS]: 


U(T;i = 0 ) = V’o(r) 
\^(r ; t) =n(T-t)^(r-t). 

Time evolution is given by the Smoluchowski operator 


N 

n{r;t) = ^2 d i 

i —1 


di 


Fi (r) — u solv (ri, t) , 


where di = and we have used dimensionless units by setting = 1 and 
ksT = 1. Here, contrary to Ref. fill , [itl [lil . [~T3 | , the initial probability 
density V’o(r) need not be the equilibrium distribution: the system can be 
prepared in an arbitrary configuration. 

Equation [2] is formally solved by 


V>(T;f) 



Q(T]s)ds 


V’o(r), 


( 3 ) 
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where e+ is a time-ordered product (see Appendix A of Ref. 0 )- At time 
t, the evaluation of a many-body function g reads 

(g) t = f g(r)mt)dr (4) 

Instead of having the probability distribution i/} evolve in time, as in 
Eq. [3j a dual formulation is sometimes preferable, in which (by means of 
a partial integration of Eq. H if is kept constant and the dehniton of g 
evolves with time, analogously to the switch from a wavefunction-evolving 
Schrodinger representation to an operator-evolving Heisenberg representa¬ 
tion in Quantum Mechanics, viz ., 


U(r;f) = V’o(r) 

\d t g(T-,t) =rt(r-t)g(T-t), 

where fit(T;i) = YliLi [& + CO + t> solv (rj,t)] • d\ is the adjoint of the 
Smoluchowski operator, with the formal solution 


/-p /o nt ( r ; s ) ds /p n \ 

g{T;t)=e_? g(r;0), 


( 6 ) 


where e_ denotes the negatively ordered exponential [12j. 


2.2. Auxiliary frame and recovery of an advection term 

Microscopic observables depend on space, via their point of evaluation r : 
g(T m ,t) —> g(v,T;t). But the prescribed velocity held generally differs from 
zero at r, so that the evolution of g(r) mingles an intrinsic evolution of the 
system and an advection by the how held. In previous studies, for instance, 
Ref. 0 , the consideration of a strictly homogeneous system (with vanishing 
spatial gradients) rendered a disentanglement of the two effects unnecessary 
and no advection term appeared in the equations. Yet, in the presence of any 
heterogeneity, such term is expected on physical grounds and is crucial for 
the study of perturbations, hence, instabilities. Here, we purport to carefully 
establish its recovery. 

To disentangle advection and intrinsic dynamics, it is helpful to observe 
the dynamics in a frame that moves with the solvent velocity at the point r a 
and time t a that will be of interest. Thus, we introduce new, time-dependent 
coordinates 

r'[r,t\=r - (r Q (t)-r Q ), (7) 

with the backward transform 
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r[r',t\ = r' + {r 0 {t) - r a ), 

where r Q (t) is the pathline of the (non-singular) solvent velocity field v solv (r, t) 
that ends at r 0 at time t Q , i.e., 


d t r 0 (t) = v solv (r D (t),t) 
r 0 (t 0 ) = r 0 . 


(In Appendix A we propose an equivalent, alternative approach, rooted in 
operator formalism rather than change of frame). At a fixed point r' in the 
new frame, g evolves with time as follows: 


D t g (r [r',t] 


Y-f) = lim 9 f + dt \ + dt ) ~ 9 ( r W, *] i r ; *) / q n 

= ^ f (r ;t)g (r [r',t] ,T;t) + d t r Q (t ) • d r g ( r [r',t\ ,T;i) 

= fl f (r -t)g (r [r',t ] ,T;t) + v solv (r 0 (t), t) ■ d r g ( r [r',t\ ,T;t) . 


Next, we notice that commonly used observables, such as the stress or the 
density, do not depend intrinsically on space, i.e., there exists a function g 
such that g(r,T) = g{r,r \,..., rqv) = g(r i — r ,..., rjv — r). Consequently, 


d r g(r,T) = d r g(r! - r,..., r N - r) (10) 

= ~X]^( r i ~r, ■ ■ ■ ,r N - r) 

i 

= -J2 di9 ( r ’ r ^ ( n ) 

i 


Inserting this result into Eq. [9l we get 


D t g (r [r', t] , T; t) = [e£i [$ + F< (r) + u solv (r i? t)] • d- 

_^s° lv (r 0 (t),t) • Ei dig (r [r', t] , E; t) 

= E£i [di + Fi (r) + (u solv (n,t) - u solv (r 0 (t),t))] • d, 


g (r [r', t] , T; t) 
9 ( r [r',t] ,T;t) 


Denoting by a prime the functions expressed in the new frame, i.e., 
f'(r'[r,t],t) = f (r,t) for a generic function f, and remarking that co¬ 
ordinates in the original and new frame are in a one-to-one correspondence, 
we arrive at 
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( 12 ) 


D t g (r [r',t] ,T;t) = d t g' (r' ,T' ;t) 

= nt'(r,ty (r'.rji), 


for any r' in the domain, where 

N 

t) = J2 [d'i + Fi (r') + Vf in', t)]-d[ 

i= 1 

and 

v' (r', t) = v soW ' (r', t) - v soW '(r 0 , t ). (13) 

Thus, an observable g' evaluated at fixed position in the auxiliary frame 
displays dynamics identical to those of its counterpart g in the original frame 
(see Eq. 0, except that the velocity field •u solv entering the Smoluchowski 
operator for g is replaced by a new field v' for g', which vanishes at r 0 . Using 
Eq. [13j we see that the evolutions in the two frames are related by 


d t g' (r',r';i) 


r'=r'[r,t\ 


= d t g (r,T;t) + v solv (r a (t),t) ■ d r g (r,T;t). (14) 


What is the advantage of switching to the coordinates in the auxiliary, 
then? First, for any (potentially time-dependent) evaluation point r' , the 
new Smoluchowski operator O' '(T 7 , i) is insensitive to global, potentially 
time-dependent translations in the original frame, i.e., offsets of the veloc¬ 
ity field v solv (ri,t\). Accordingly, it only depends on the velocity gradient 
Ka/3 ( r ,t) = d/ 3 V^ lv ( r,t ) = d'pv' a (r'[r,t],t), responsible for the deformation 
of the structure. A second considerable benefit emerges for the specific point 
of evaluation r' = r 0 , where the field v' vanishes at all times. (Note that 
in the original frame this point moves as time passes). Then, in Eq. [bU the 
effects of advection and of the intrinsic dynamics are clearly separated: the 
latter are reflected by the evolution of g' at the point r' = r a in the auxiliary 
frame, where there is no flow, while the second term on the right-hand side 
(rhs) is the desired advection term; let us once again emphasise that this 
non-local term is physically crucial in a heterogeneous system. 


2.3. Leading-order locally homogeneous flow 

Formally, Eq. [6] conveys the impression that a microscopic observable g , 
albeit evaluated at a given point r Q and time t a , depends on the configuration 
T of all particles throughout space , and not only at r = r a , and hence requires 
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the knowledge of the whole solvent velocity held -u solv (r, £ ^ t Q ). However, 
making use of the short range of usual observables, we purport to bolster 
the intuition that, to leading order, (g(r o))^ o is mainly determined by the 
history of the velocity gradient ft(ro (£),£) along the solvent pathline r 0 (£), 
with r(t a ) = r a . 

We dehne the range of a microscopic observable g{r) as the distance 
beyond which the particle configuration becomes irrelevant. More precisely, 
range (g) is the minimal radius of a disk T> centred at r such that, for any 
two particle configurations r^ 4 ) and F^b coinciding over T >, he., such that 

r ?: (A ) = if G V or r^ G V, (15) 

g{r,T^) and g(r,T^) are equal, to a good approximation. For example, 
the range of the density observable p(r) = <5 ( r ~ r j) is 0 + and that of 

the stress er(r) is bounded by the cut-off distance of interparticle interactions. 

If the range of an observable g{r Q ) is small compared to the lengthscale 
l rvj k)\7k over which the velocity gradient varies, we are tempted to replace 
the global inhomogeneous flow with a much more tractable affine (he., ho¬ 
mogeneous) velocity field that coincides with the inhomogeneous one around 
TV 

This comes down to approximating the genuine Smoluchowski operator 
(appearing, e.g., in Eq. [6]) with 


a 


N 

Lm( r ;0 = ^q( r ) +5Z [ vSOlV ( r o,t) +K(r a ,t) ■ {Vi - r a ) 


%— 1 


di 


where 


N 


^ q (r)^E^ + F *( r )]-ft. (16) 

2=1 

How large is the error due to this approximation? At time £, the error reads 

— el u .. I5f(r 0 ,r). 


Jo^(s)ds Jo^Lm( S ) ds 


In particular, at £ = 0, it is zero, and the first-order term in £ yields its initial 
growth rate, 


rt ( N 

fit ( s ) “ °Lm( s ) g(r 0 ,T)ds = O PT I|Vk|| \\n -r* 0 || 2 \dig{r 0 ,T)\ 


'0 L 


^ 2=1 


= O ( A r x>(r) ||Vk|| |range(g )| 2 max \\dig(r 0 , T))|| 

ie{i,...,iV} 





where Nd(T) is the number of particles within disk T> in configuration T. 
Clearly, ||V/c|| |range(< 7 )| 1 2 arises because of the local deviations from affinity. 
The second-order term (quadratic in fit) in the expansion of the approxi¬ 
mation error also contains contributions in ||V/c|| |range(g)| 2 ; some are mul¬ 
tiplied by fijq (which tends to restore the equilibrium configuration), while 
the others involve ?; solv (r Q ,s), which drags particles away. Indeed, through 
advection with the solvent velocity, some particles, which initially lay far 
from r Q (where the affine approximation is poor), will enter the region T> 
where they become relevant for the computation of g. Consequently, for the 
approximation to work best, particles close to r a should move as little as 
possible. This is exactly why it is advantageous to switch to the auxiliary 
frame introduced in the previous section: the auxiliary driving field v' ( r',t ) 
vanishes at point r' = r a at all time&Q. In that frame, the approximate 
evolution is ruled by 


with 


</(r 0 ,r';fo) = e_ 


_ Jo° n L( r '; s ) (is 


g(r a , r' ; o), 


(17) 


N 

^ m (r , ) t) = ^2 [d[ + Fi (T') + k' ( r 0 ,t ) • (r/ - r 0 )] • d[. 

1=1 

Recalling that, by definition, the original frame and the auxiliary one coin¬ 
cide at time to and thus (g'{f‘ 0 )) to = (g{r 0 )) to , we can now come back to 
the original frame, using the change of coordinates of Eq. [71 and confirm the 
intuition that, to leading order, (g(r 0 )) t is governed by the velocity gradi¬ 
ent along the pathline, i.e., {k (r D (t), t)}. This is the intrinsic part of the 
dynamics. As a reminder, the time derivative of ( g(r 0 )) t also involves an 
extrinsic part, namely, the advection term in Eq. [141 

To proceed, physically motivated approximations, expressed as projec¬ 
tions onto relevant variables, are performed onto the intrinsic dynamics. In 
the end, these approximations, conducted in Fourier space, shall heavily rely 
on the possibility to treat the driving flow as (almost) globally homogeneous, 
whereas the system under study may be globally very heterogeneous. The 
problem is solved by performing these approximations in the homogeneous 
auxiliary system of Eq. El which is a reasonable surrogate for the origi¬ 
nal one if the flow is locally homogeneous. This is a first step towards a 


1 Note, however, that the quality of this locally homogeneous approximation will dra¬ 
matically worsen with the duration of the memory of the system and the magnitude of 

K. 
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systematic expansion of the velocity field in the auxiliary frame, starting 
(in this paper) with a uniform velocity gradient n(r,t ) = n(ro(t),t), then 
considering the gradient Vk of k at r q (t), etc. 


3. Projection Scheme 

3.1. Sets of slow variables 

In a typical mode-coupling spirit, the slow intrinsic evolution (with re¬ 
spect to microscopic time scales) of a generic observable g such as the density 
or the stress will be captured via its projection onto [i.e., cross-correlation 
with) familiar slow modes. The observables will be expressed in Fourier 
space, where the collective dynamics are best captured. Since the global 
density p q =o (where q represents a wavevector in Fourier space) is the only 
conserved quantity in the problem, i.e., dtp q =o = 0, and the relaxation time 
of p q diverges in the limit of small q, we define the linear density modes 
|p<j) Q £ R rf } hr Fourier space as a first set of slow modes, associated to the 
projector Pi, viz., 

A = ^~]Pq) atq ( Pq ’ 
q 

where S q = IV -1 {p q p q ) is the static structure factor, and its complemen¬ 
tary part Q\ = 1 — P\. It should be noted that the ensemble average in 
the projection is performed with respect to the equilibrium distribution ip eq 
(denoted by (•) here), whereas averages over the initial distribution ipQ shall 
be denoted by (-) 0 , 


A5(T) 



p q (r') g (r') ■i/’eq (r') dV 


In Ref. jl2| . Brader et al. noticed the absence of any coupling with linear 
density modes in a purely homogeneous flow and thus further projected the 
dynamics onto density pairs {p k p q i k, q € R^} with the projector 


P ‘2 = J2 

k> q 


1 


N 2 S k S, 


k^q 


' (, PkP* 


q> 


where the Gaussian approximation (p k P q PkPq) ~ {p k Pk) ( P q Pq) = pr2 ‘S k S q 
was used, and its complementary part <52- Although this section comes in 
the wake of Ref. [12l |. we shall not neglect the couplings with linear density 
modes from the outset, because the flow is not strictly homogeneous. 
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3.2. Generalised Green-Kubo relation 

To prepare the projection, we recast Eq. [ 6 ] into a form which better 
highlights the deviations from the initial configuration occurring throughout 
the past. 

In the Schrodinger-like formulation, we denote by Sip these flow-induced 
perturbations, viz., 

ip(T,t) = ip 0 {T) + Sip(F,t). 

Since dtip(T,t) = 17(T,7) [V’o(r) + <5^(T,t)], solving for Sip yields 

G f! ntr,s)ds 

ip(T,t) = ipo(T) + / dt ie + 41 n(r,ii)V-o(r). (18) 

■Jo 

For the time being, 17 is the Smoluchowski operator for a generic flow field 
u solv , but the approximations performed in the following (see Section ED 
will hinge on its being close to homogeneous, so one might already think of 
17 and 1ft as their homogeneous auxiliary-frame surrogates of Eq. [TT] that 
is to say, mentally consider the evolution in the auxiliary frame, with the 
replacement 17 1 —> J7^ . 

Applying Eq. [18] to an arbitrary observable g, e.g., g = a, and partially 
integrating the time-ordered exponential, we arrive at a generalised Green- 
Kubo (gGK) relation, expressed in the Heisenberg-like representation, 


(g) t = / dTg (T) 


f G f* Q(r ,s)ds 

^o(r) + J dr j^ dhe^ 17 (r,t!)^o(r) 


P / , (J Clps)ds 

{g) o + J dti /l7 t (7i)e_ 1 g 


(19) 


The Green-Kubo nature of Eq. [T9] becomes clearer if the integrand is 
rewritten as 

f [J nt ts)ds 

J g. 

n(T,ti)ip 0 (T) is thus the deviation from ipo(T) created at time t\ (per unit 
time). For instance, for simple shear flow, starting with ipQ = ip e q , the 
deviation couples strain rate and stress: 17(T, ti)ip eq (T) = 7(7i)cr E < ; (r)'i/7eq(r), 
where a xy is the shear element of the Kirkwood stress tensor and 7 (t) is the 
imposed shear rate. 


3.3. Projected dynamics 

Let TJJ(t,ti) = eJ f -i dt2fl ' t2> be the propagator appearing in gGK (Eq. [T9]), 
associated with the full dynamics 17^. We split U t into a part u\(t,t\) = 
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Figure 1: Schematic diagram of the decomposition of the full propagator 
W(t,t\), associated with the operator fit, as performed in Eg. 1201 


g/tj dtxQiCi (*i) evo l ves purely orthogonally to Pi and a part that inter¬ 
acts at least once with P\ (the notation t 2 referring to the time of the last 
interaction, see Fig. [[]), viz. 

u\t,h)= f dt 2 u\t2,t 1 )p 1 n\t2)ul(t,t 2 ) + ul{t,t 1 ). (20) 

Jt\ 

Inserting the decomposition of Eq. [20] into gGK (Eq. [T9l) . we arrive at: 


{g) t - (a )o = 


[ dt 2 [ dti (^(ti)U^(t 2 ,t 1 )Pin\t 2 )ul{t,t 2 )g 

Jo Jo ' 

+ [ dti (&(ti)ul(t,ti)g\ 

J o' ' 0 

ft rt2 , v (p*Stf(t2)ul(t,t 2 )g 

/ dt ^y2 / dti/n\ti)u\t2,ti)p q 

Jo </ h _1_ 


NS„ 


+ 


( Pl)t2 (pl)o 

J dti (n\ti)ul(t,ti)g^ 


( 21 ) 


(Orthi) 


where we have made use of gGK (applied to density fluctuations, g —> p q ) in 
the last equality to reduce the first brace to (pq) t — (Pq) 0 - One thus arrives 
at: 


{g) t - (a)o = 


/<#>£ 


Pq) 


t2 



(p* q ^(t 2 )ul(t } t 2 )g^ 

NS~q 


+ (Orthi). 


( 22 ) 
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First, we focus on the dynamical correlator {t 2 )tj\(t, 1 2 ).9y on the 

rhs and introduce the identity P\ + Q\ = 1 as follows: 

'p*ni(t 2 )ul(t,t 2 )g / ± , , 

(t 2 )ul(t,t 2 ) (Pi + Qi)g 


NS n 


NS q 


'E v £ M St & 1 *) 


+ 


NS q \ rq 


P*q^(t 2 )ul(t,t 2 )Qig\ , 


where the vertex V£ = quantifies the coupling of the observable g to 

the density mode p in the equilibrium distribution and 


"lid M2) = 


NS„ 


P*q&(t 2 )Ul(t,t 2 )p k 

NS~„ 




(23) 


+iq ■ ( -Fq + Vq*(t 2 )) u\ (t, t 2 )p k 


is a memory kernel evaluated in the equilibrium distribution, with 


F q = 


Y F je iqrj and v q (t) = Y V 

3 


solv. 


rj,t)e 


-iqrj 


(24) 


3-4■ Application to the density observable 

Before turning to our main interest, he., the stress, we wish to illustrate 
the principle of the projection scheme for a generic flow, but on a simpler 
observable, namely, the density g = p p , p £ M d , for which the complement 
Qig vanishes by definition. The following calculations need not be performed 
in the homogeneous auxiliary frame; they hold true for an inhomogeneous 
flow. 

Applying Eq. [22] to density modes (V£ p = = bk, P ) leads to 


(Pp)t (Pp) 0 — 


[ dt 1 Y [Wti “ Wo M <JP (*»*i) + (Orthl), 

3 0 r. 


where M qp (t,ti) is given in Ea.l23land (Orthl) = f* dt\ (ti)ul(t, ti)p p \ . 
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Taking a derivative with respect to time t yields 
dt (Pp) t = — [(Pg)t ~ (Pq)o] MqJ (t,t) 

Q 


■f*i E 

Jo „ 


*1 


(Pq) c 


d t Mqp (t, ti) + 5p(i). 


(25) 


Here, we have used the explicit notation S p (t) for <9j(Orthl); 

rt 

Sp(t) =-ip ■ (F p + v p - ipp p ) o- / (-F’p + 'Up))o- 

Jo 

(26) 

The term Mgp (f,f) = — (IVS'g) -1 (p*f2^(t)p p ) can be simplified. Using 
the equilibrium (■ i.e p solv = 0 ) Smoluchowski operator flj q from Eq. [TO] we 
can write 


(p*H f (t)pp) = ^ p*qtti q pp ) 

where, using partial integration, 




The second term, 





q / / 



-p 2 5 qp N. 


(\Pq l,solv ’ djPp^j = -Pgp? is most easily simplified 
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by first backward-Fourier transforming it with respect to p only, viz., 

f q ( r o) = Y vSOlv OicO • d iP (n>, 

= ~(pqY1 wSOlV ( r J>*) • d ro s (ro - rj)j . 

= -9r 0 - (pqYI vS ° W ’ *) 5 ( r ° _ ^ 

= - (p^ solv (ro, i) • <9 ro p (r 0 )^ , 

where we have used the incompressibility of the velocity field, viz., d ro ■ 
■ysolv (ro,t) = 0. It suffices to transform F q (ro) back into reciprocal space 
(with respect to r q) to obtain: 

F qp = — 'y ^ ( Pq v k{t ) ■ i (p — k) Pp—k ) 

k 

= -v k (t)-iq{p*p q ) 

= -iNv s °}_ v q (t) ■ qS q . 

Collecting these contributions into Eq. [25j one arrives at an equation of 
evolution of density fluctuations: 


dt {pp) t + Y v p-q (*) ' iq (P*)t = < Pp)t ~ Y [ dt i d t M qp (t, h) ( Pq ) t 

„ &P „ Jo 


+S p (t). 


(27) 


Let us emphasise the physical meaning of the different terms: 

(i) the second term on the lhs is the advection term established in Sec¬ 


tion 

(ii) the first one on the rhs permits the relaxation of fluctuations through 
diffusion. Note that the normalising factor is expected, because the 
relaxation of a density mode p p does not require single-particle diffusion 
over a lengthscale ||p|| -1 ! Besides, if the density field is smoothed via coarse- 
graining (i.e., p p —> p p (f> p , with (j) p akin to a Gaussian of half-width a few 
particle diameters), the fast relaxational modes at high p are suppressed, 
and the normalising factor then tends to Sq 1 , which is directly related to 
the compressibility of the suspension. 
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(iii) The second term on the rhs reflects the evolution in orthogonal space 
of the density fluctuations created in the past and their final coupling back 
to present fluctuations. 

(iv) S p (t ) is a source term that results from interactions with nonlin¬ 

ear density modes. Consistently with the expectation that density hetero¬ 
geneities in an incompressible flow are due to collective effects, e.g., stress 
equilibration, it is the only term that can potentially create density inhomo¬ 
geneities, insofar as the other terms are associated to pre-existing density het¬ 
erogeneities. If a homogeneous flow is imposed to an initally uniform system, 
translational invariance imposes that, for finite p, (^Fp^ = 0, {p P ) 0 = 0, 

ip ■ (v p ) 0 = 0 in Eq. [26] ergo S p (t ) = 0 at all times, for all p ; as expected, 
the source term then vanishes. 

Finally, by comparing Eq. [27] to the mass conservation equation, 
dt{p P )t + ip- (jp U )t = 0 , 

where j™ 11 = \ rjer iv ' r i is the colloidal flux, it can be seen that, in a 

heterogeneous flow, the colloidal velocity u coU (r,t ) = generally 

differ from the driving solvent velocity D solv (r,t). 

3.5. Orthogonal dynamics 

Let us come back to a generic observable g and refine the description 
of term (Orthl) in Eq. [2TJ In principle, the propagator decomposition of 
Eq. [20]can be iterated, and the propagator u\, split into a part evolving with 
P 2 Q\0) and an orthogonal part and so on, ad libitum. Schematically, 
one would then get 

(d)t ~ (9)0 - [ dt 2 ( Wt 2 _ Wo) (•) M qk (*» *2) ( 28 ) 

~ f dt 2 (( PkPp) t2 ~ (PkPp) 0 ) ^2 (') t 2 ) + ..., 

0 k,p,k',p ' 

and the orthogonal evolutions denoted by M^ n \t,t 2 ) would then be con¬ 
strained to smaller and smaller spaces. However, following Ref. ji2l |. we 
adopt a more pragmatic approach by directly introducing the projector P 2 
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in (Orthl), viz. 


(Orthi) 


where 



dti (n\ti)QiUl(t,ti)g^ 

dt 1 (n\t 1 )Q 1 P 2 ul(t,t 1 )P 2 g) Q 



(ti)^ 


PiPpUi{t,h)p k >p pl 

N 2 S k S p 


k'.p' 


k f >p f 


( 29 ) 


^0 ,k,p (^l) 



(t i) Q i pk Pp ^ ^ 

(pI'P’p's) 

ms k .s p , 


are vertices that represent, respectively, the creation and relaxation of bilin¬ 
ear density modes with respect to the initial configuration and the coupling 
strength of the observable g to these modes. If the initial configuration ipo 
and the velocity gradient at t\ are close to homogeneous, then Vq^ p (ti) is 
nonzero only for p ~ — k, and if g is a spatial average, that is, if only small- 
wavenumber modes contribute to its Fourier decomposition, then V£, p , is 
nonzero only if p' ~ — k'. 

To conclude, the density-pair correlation function \P k p p Ul(t,ti)p k /p p /\ 


is approximated through Gaussian factoring, which is a central approxima¬ 
tion of mode-coupling theory [21|: 


{plPpUt{tM) Pk 'P P >) ~ (plp* p u\tM)p k 'P P ') 

N 2 S k S p ~ N^S k S p 

(pl ui {t,h)p k f) (pp£/ f (Mi)/y) 

NS k NS P 


where we have introduced the transient density correlator 


$ qk{ti ^ 2 ) 


(, PqU\t,t 2 )pk) 
NS q 
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This correlator indicates how fast density fluctuations relax, in the presence 
of flow. A central result of the rheological extension of MCT by Brader 
and colleagues 12] ascertains that, for a translationally invariant (i.e., ho¬ 
mogeneous) flow, & q k(t,t 2 ) is non-zero only if q coincides with k in the 
solvent flow frame, i.e., if the flow advects wavevector k at time t <2 into 
J t * dsn(s) 

q = k(t,t 2 ) = e + 2 at a later time t, where n is the uniform velocity 

gradient (recall K a p ( r,t ) = <9^r;® olv ( r,t )). Interestingly, the norm of k(t,t 2 ) 
increases with the deformation; consequently, <f> is effectively evaluated at a 
smaller lengthscale, where thermal relaxation occurs faster. But flow het¬ 
erogeneities induce additional cross-couplings between p* and U^(t,t2)pki 
through the interaction between the structure (the density modes p q and 
Pk) and the flow field (the velocity gradient modes /%>(s))cl 

After collecting all terms, we arrive at the following expression for Eq. [2TJ 


ii (Pq) o 
1 


<S>t - (s)o ~ ~ [ dt ±Y^ 

J o q 

+ NS W, 


L k 

t 




+ 


I 1 1) Iqy 7 ^l) ^ pp' (b ) ■ t ■ 1 ^) 


k>p 

k r >p f 


We should pay attention to a possible issue with Eq. [30} if one is not 
cautious, there is a risk that the uncontrolled approximation of (Orthl) may 
create spurious inhomogeneities that violate translational invariance, even 
in cases where it should theoretically be obeyed. We shall ward off this risk 
in a pragmatic way by making a judicious choice for the transient density 
correlator and ensuring the respect of fundamental physical principles in the 
final equations. 


2 This can be seen by replacing U*(t,t 2 )pk in & q k{t,t 2 ) with its gGK expression 
(Eq. 1191) . The resulting integral involves a rate of deviation fl(s)ift eq = — Ej vSoW ( r j,s) ■ 
Fjip eq , which is proportional to Ko(s) : <xo if the flow is homogeneous, but comprises 
non-zero Fourier modes otherwise. 
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4. Severe approximations lead to a constitutive equation of the 
White-Metzner class 

Armed with the projection scheme of the previous part, we are now 
capable of studying the stress observable cr q . Direct application of Eq. 1501 
gives 


( <J q)t (^q ) o 



dt i 

Q 



X 


(Ml) + jA- 



(*:i) V k' Q p '®kk' (Ml) <$> pp > (t, 11)■ 


(31) 


fc / >p / 


The usual MCT protocol involves the derivation of equations for the 
transient density correlators $> qp (t,ti), with the help of the Zwanzig-Mori 
projection formalism. However, in the presence of flow heterogeneities and 
in non-trivial geometry, this would be a difficult task, which we bypass here: 
bearing in mind our main goal, namely, a study of the visco-elastic insta¬ 
bility, we opt for a drastic simplification of the equations. We shall thus 
perform very strong, and mostly uncontrolled (but explicitly mentioned) ap¬ 
proximations and we do not expect them to preserve the quantitative details 
of the full theory. Nevertheless, we aim to arrive at a tractable model that 
displays shear-thinning and correctly captures the low-shear-rate and high- 
shear-rate regime of the flow. This will allow us to test the phenomenology 
of our extension of the framework to inhomogeneous situations and flow 
instabilities. 


4-1- Reduction to a locally homogeneous flow 

To start with, we assume that the flow only moderately deviates from 
homogeneity, so that the locally homogeneous flow approximation of Sec¬ 
tion 12.31 is valid. Under that assumption, the source term S p (t ) for the 
density in Eq. [27 ] vanishes, and no density inhomogeneity is created (effects 


such as shear-concentration coupling [22l. 12 31] are thereby precluded). Thus, 
we suppose (p q ) t ~ {p q ) Q ~ N5 q p. 

Furthermore, if the velocity gradient varies little over the distance trav¬ 
elled by a material volume before its microstructure relaxes, the intrinsic 
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dynamics of the local stress <r(r,t) are governed by the history of the lo¬ 
cal velocity gradient {k t}. In other words, all non-local effects 

due to, e.g., stress waves emitted by distant regions are discarded. Within 
mesoscopic regions of the material, translational invariance is obeyed, so 
that only one cross-coupling subsists in the transient density correlator, viz ., 

Under these assumptions, the vertices in 

Eq. [2U boil down to [12j 


v«, P ' = fcV 


or 

'V 


N 2 k'Su/S. 


k' D p' 

S' 


Jk f ,—p f 5 


Thus, we recover the formula for homogeneous flow derived by Brader et 
al. |15|, |12|], namely, 


a(t) = - 


[ dt\ V A(k',t,ti) 

Jo 


-fi- i k ' ■ ■ k') 

(32) 

ft dSK,(s) 

where where we have introduced the Finger tensor B{t,ti) = e, 1 
dsK T (s) S', S', 

e_ x and A(k,t,ti) oc ^ ’ V collects all relevant equilibrium 

K 

properties of the material. 


4-2. Schematic approximation 

To proceed, we follow the schematic approximation conducted in Ref. 
by dropping all existing wavevector dependences (or, better said, by focusing 
on the most relevant wavevector) in Eq. | 


15] 


viz., 


cr(t) = -v a [ dt\ —— 77 r ~ *£ 2 (UU), (33) 

Jo Ot 1 

where v a ss A^Tn (with n the number density) sets the scale of stress fluc¬ 
tuations. Partial integration yields 

a(t) = v a fdhB^t (34) 
Jo ot] 
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4-3. Approximation of the transient density correlator 

At rest, the thermal relaxation of density fluctuations takes a time r a 
that, according to MCT, diverges in the ideal glassy phase. But the pres¬ 
ence of a solvent flow distorts the material structure and accelerates this 
relaxation. Arguably our crudest assumption will now consist in proposing a 
“phenomenological” characteristic time scale r for the long-run decay of the 
transient density correlator, which takes this flow-induced relaxation into 
account: 


d . . 

M* itM = 

= 

with r [e] = 


$(Mi) 

2 T [e (i)] 

$(Mi) 

2 t [e (ii)] 

_ fa _ 

1 + 2 ar a \JJ 2 (e)’ 


(35) 

(36) 

(37) 


Here, e (t) = K, C)+^ ^ jg the strain rate tensor, J 2 is the second tensorial 
invariant of a deviatoric tensor, viz.. J 2 (e) = (this reduces to J 2 (k) = 

|y 2 under simple shear), and a is a material parameter that quantifies shear¬ 
thinning. More precisely, a~ 1 is an “inverse yield strain” that describes how 
much strain is required to erase the memory of the local structure fldl j . 
The following limit cases are enlightening: at vanishing shear rate, r (e) 
tends to the quiescent relaxation time r a , while at high shear rates one gets 
r (e) ~ (ccy)" 1 (under simple shear). 

Equations I35ll36l lead to 


<f> 2 (t, 1 1 ) = exp 



(38) 


4-4- Constitutive equation 

With these severe approximations in hand, we now differentiate Eq. [M] 
with respect to time t. 

dt t [e (t)] r [e (t)] 

where I is the identity matrix and we have used the equality = 

ti=t 

r[e(i)] _1 as well as 9B gf 1 ^ = re(i) • B(t,t \) + ■ n(t) T . Splitting 
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a into a quiescent (ie., n = 0) pressure part v a I and a driven part <x d in 
Eq.Ei we arrive at 


d 

dt l 


(t) = V a K ( t ) + K T ( t ) + K)(t) • CT a {t) + <T d (t) • K,(t ) 


T 


g d (fl 

r [e (t)]' 


(40) 


Finally, one should recall that Eq. 00] deals with the intrinsic dynam¬ 
ics, i.e., the evolution in the (homogeneous) auxiliary frame introduced in 
Eq. 03 The full evolution of the observable in the lab frame is recovered by 
applying Eq. 04] which restores the inhomogeneity advection term v ■ Ver 
found in Section 12.21 responsible for the advection of stress fluctuation with 
the microstructure, hence with the flow. After dropping the “d” superscripts, 
we thus obtain 


2) , , <x(r,f) .. , 

^ (r ' t)+ 7yM)r 2ive(r ' t) ' 


(41) 


where the upper convected derivative (a.k.a. advected Truesdell derivative) 
is defined by 


2)<t da „ . . . . , . , , T 

Remember that the final relaxation time depends on the intrinsic dynamics 
and the shear rate and is given in Eq. [37] 

Interestingly, this relatively simple constitutive equation belongs to a 
class of (mostly) phenomenological models initially put forward in the con¬ 
text of polymer melt rheology by White and Metzner (WM) 0 , on the basis 
of symmetry considerations. It obeys the principles of locality, causality and 
material objectivity. 

In addition to the configuration-based stress a, which is associated to the 
colloidal microstructure, we include a Newtonian contribution to the stress 
accounting for viscous dissipation; the latter can indeed become significant 
at large shear rates. The total stress then reads 


£ = <t + 2 rj s 


(42) 


Note that this model has already been studied by Papenkort and Voigtmann 
in the case of a channel flow 
played no role. 


24| , where however the advected derivative 
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4-5. Bulk rheology and parameter fitting 

Within the WM-type model, we first consider the rheological properties 
of the bulk shear flow, he., k = 7 *ei , prior to any potential instability. 
From the constitutive equation, using tjq as a shorthand for v a T ai it is easy 
to show that 


Oil 

_★ 

°12 

°22 


0 


mi 


1 + ar a \i 

2?7oT( 


7* 


1 + crr Q |t* 


(43) 


Turning to linear rheology, the storage and loss moduli associated to the 
WM-type model are identical to those of a Maxwell model with a Newtonian 
contribution, viz ., 


G'(u) 

G» 


VO T a LO“ 

1 + W 2 T 2 

Vqu 

1 + w 2 r 2 


+ VsU. 


To fit the model parameters v a , r Q , ct, and ? 7 S , we compare the ct* 2 ( 7 *)- 
flow curve and, with less emphasis, the storage and loss moduli to experimen¬ 
tal measurements by Siebenbtirger et al. on suspensions of ~ lOOnm-large 
colloids with a thermosensitive PNIPAM shell (which affords a sensitive con¬ 
trol of the effective volume fraction through the tuning of the temperature) 


25]. In Ref. [25], a schematic version of the MCT equations was shown to 


provide excellent agreement with both the oscillatory shear and the steady 
shear measurements, for several effective volume fractions across the glass 
transition, while the solutions of the microscopic MCT equations were tested 
in Ref. [lj] and yielded a consistent first principles description. 

Figure [2] presents the flow curve fits obtained with our much cruder (but 
also much more tractable) constitutive equation; there is no doubt that the 
quality of the fits has suffered from our strong approximations: this con¬ 
firms that the memory effects encoded in <f>(i,fi) are more subtle than our 
simple Ansatz of Eq. [37] Nevertheless, the downward bending of the stress 
at low shear rates, which originates from thermal relaxation, and the strong 
shear-thinning effects, as well as the viscous behaviour at high shear rates are 
reasonably well captured. Let us also mention that our fitted model param¬ 
eters (see Table [T]) have values comparable to the corresponding quantities 
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4>efi 

T a 

a 

Vo 

Vs 

0.519 

1.7 

22 

9 

.0 

1.2 

0.600 

900 

35 

8.1 

10 a 

1.3 

0.626 

6.0 • 10 4 

30 

5.4 

10 b 

1.5 

0.639 

8.0 • 10 5 

30 

6.4 

10 e 

1.9 

0.641 

3.0 ■ 10 b 

40 

4.8 

10 v 

2.0 


Table 1: Model parameters used to fit the experimental flow curves in Fig. [2] 
and linear spectra in Fig. [3]. 

^3 7 m 

As in Fig. [2] stresses are in units of and times in units of B . 


in Ref. [25| : note, for instance, that t q increases dramatically with 0 e ff and 


that a is roughly of the order of the inverse yield strain. 

We bestowed less importance to the linear rheology moduli G' and G" 
in fitting the model parameters, because the rest of the paper deals with 
steady shear; still, Fig. [3] shows that the experimentally observed trends are 
also present in the model, although there is no quantitative agreement. The 
origin of the deviations is the broad distribution of relaxation times neglected 
in the White Metzner model but contained in the MCT description. 


4-6. Base flow in Taylor-Couette geometry 

Now, we focus on the base flow in a curved geometry, and more specif¬ 
ically in the Taylor-Couette cell used by Siebenbiirger and colleagues in 
Ref. (25], In such a rheometer, the fluid flows in the annular region be¬ 
tween two co-axial cylinders, of radii r* = 13.33 mm and r D = 14.46 mm 
(relative gap width e = 0.085), due to the rotation of the inner one. This 
geometry will be kept fixed for the rest of the paper. 

Turning to the determination of the flow, we note that, unlike in Sec¬ 
tion [2] the velocity profile v(r) is no longer prescribed. Therefore, we need 
to close the equations by complementing the WM-type constitutive equation, 
(Eq. Idlll with the inertialess momentum conservation equation, 

0 = V • £ - Vp, (44) 
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(a) < p eB = 0.519 (b) <t> eS = 0.600 




(c) q ieff = 0.626 (d) <j> e s = 0.639 



(e) 0 eff = 0.641 


Figure 2: Flow curves for different effective volume fractions. (Red lines with 
crosses ) experimental measurements from Ref. (25| : (dashed green lines ) fits 
with the WM-type model (see text); parameters are listed in Table I. 

R 3 

Stresses are given in units of where Rh is the hydrodynamic radius 

of the colloids, and shear rates are non-dimensionalised with (>nrilial ^ rt r witfi 
’ k B T ’ 

rjsolv the solvent viscosity. 
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(a) 0eff = 0.519 



(c) 0 e ff = 0.626 



(b) </> e ff = 0.600 



(d) 0eff = 0.639 



(e) ^eff = 0.641 


Figure 3: Linear rheology measurements for different effective volume frac¬ 
tions. (Blue squares ) experimental storage moduli G'(u) from Ref. [iij ]; 
(dashed blue line ) fit with the WM-type model (see text). (Green dots ) 
experimental loss moduli G"( w); (dashed green line ) fit with the WM-type 
model. The fitting parameters, as well as the stress and frequency units, are 
identical to those used in Fig. [2j 
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where p is the pressure, and with the postulate of incompressibility, 


0 = V ■ v. (45) 

In the considered geometry, the base flow is purely azimuthal and has no 
dependence on 9 or z, in cylindrical coordinates. It follows from Eg. 1441 that 
the total shear stress must satisfy 

E(r) = E(rOj|. 

Inversion of Eas. 1421 and 1431 yields, for 7 ( 77 ) > 0, 


7 *( r ) = 


2 ap s 


a profile of which is plotted in Fig.[4]for a particular applied shear rate. Since 
7 *(r) = Vq —we can solve numerically for the velocity profile vo(r). It 
is perhaps worth indicating that an analytical solution exists for r] s = 0 : 


vt(r) 


ar a V v a oi - £(r) 


(46) 


Incidentally, note that the quality of the approximation p s = 0 is not fixed by 
the inequality p s <C 770 , but by the more stringent condition p s 


5. Shear-Thinning suppresses the visco-elastic instability 

5.1. The visco-elastic instability 

For almost one century, it has been known that the inertial flow of New¬ 
tonian liquids is prone to a centrifugal instability at large Reynolds num¬ 
bers (or, more precisely, at large Taylor numbers), whereby counter-rotating 
vortices, known as “Taylor vortices”, appear and break the full cylindrical 
symmetry of the base flow 0,0- Shear-thinning fluids are also prone to 
this instability j 28| . which may even be enhanced, owing to shear-thinning, 
for polymeric solutions 0]. 

This type of instability is driven by inertia and opposed by viscosity. 
In dense colloidal suspensions, the vanishing Reynolds number precludes it. 
Yet, there exists a distinct type of instability, which does not require inertia: 
the so called visco-elastic instability was first analysed by Muller, Larson and 
Shaqfeh [3|, 0, SI] and has notably been observed in polymer solutions [3] 
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0.00035 


0.00030 

0.00025 



Figure 4: Shear rate profile of the base flow in the Taylor-Couette cell, 
with model parameters corresponding to </> e ff = 0.626, for imposed (non- 
dimensional) shear rates at the rotor (blue) 7 (77) = 10 -5 and (green) 7 (77) = 
3 • 10 " 4 . 


and (semi-dilute) worm-like micellar solutions [i|, in Taylor-Couette, cone- 
a.nd-plate and parallel-plate rheometers (for a review, see Ref. |32l|). To 


grasp the importance of this finding, one should remember that, for decades, 
scientists have measured the rheological properties of polymer solutions and 
polymer melts in such setups, under the assumption of a well-defined (base) 
flow. 

The precise mechanism driving the instability is still, to some extent, 
unsettled, but it is clear that curved streamlines and material elasticity (i.e., 
a finite structural relaxation time) are vital for the (linear) instability to 
develop. Indeed, because of the curvature in, say, a cylindrical setup, the 
normal stress £00 is coupled to the radial velocity component; this can lead 
to a positive feedback mechanism, whereby a spontaneous radial velocity 
perturbation alters £ 00 , which further amplifies the perturbation [ 3 ^]. To 
assess the effect of the material’s properties on the instability, Larson used a 
model featuring a distribution of relaxation times and an adjustable shear¬ 
thinning exponent, and showed that shear-thinning tends to stabilise the 
flow, by shifting the unstable region (in parameter space) to larger applied 


shear rates [33]. Nevertheless, the author noted that “an elastic instability 


can occur even in highly shear-thinning entangled polymer solutions”, in 
the light of his calculations and experiments. Clearly, we are interested 
in knowing whether this also holds true for highly concentrated colloidal 
suspensions: Does the dramatic shear-thinning behaviour of these materials 
allow for a visco-elastic instability within the experimentally probed range 
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of shear rates, according to the model introduced in the previous section? 


5.2. Pseudo-spectral method for linear stability analysis 

To address this question, we analyse the stability of the base flow in 
Taylor-Couette geometry (see Section 14.6|) with respect to linear perturba¬ 
tions. Although many studies have focused only on axisymnretric instabil¬ 
ities, non-axisymmetric modes were shown to be even more unstable in a 
number of cases |.'14l. I.S5l|. so we investigate three-dimensional perturbations, 
i.e., perturbations 5(f with a spatial dependence not only on r and z, but 
also on 6, viz., 


6cj>(r, 6, z,t) = ( 5a rr 5a r g 5cr rz 5agg 5crg z 5a zz 5v r dvg Sv z 5p ) 


T 


We resort to a pseudo-spectral method. 

(i) We start by linearising the equations of the problem, comprising the 
six constitutive equations (Eqs. [4T|), the three momentum conservation equa¬ 
tions (Eqs. l44|h and the incompressibility postulate (Eq. ld5D . around the base 
flow so as to obtain a system of linear equations of the form 



\ 


V 

V _ 


dt 

o 

o 

0 


■V 

A 


5<f> 


0 / 




Then we transform the endomorphisnrs A and CA into matrices, viz ., A —> A 
and C* —> L*, with the following generic procedure: 

(ii) 5cf> is Fourier-transformed both in the azimuthal direction (wavenum¬ 
ber m ) and in the axial direction (wavenumber k ), viz ., 




E 

1=-00 


5(j>(r, m, k, t ) e 


,im9 ikz 
& 1 


111 


is discretised along the radial coordinate r by sampling its values 
at the Chebyshev-Gauss-Lobatto points r n = 1.0 + |[l-|- cos (tt^)] , 0 ^ 
n < M, for a given M £ N*. We will typically use M ~ 2 6 interpolation 
points across the gap. Radial derivatives d r 5cj) are then written as matrix- 
vector products of the form D r ■ 5<f> 361 ]. 
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(iv) 5(f> is Laplace-transformed with respect to time, with Laplace coor¬ 
dinate s, viz., S4>(r,m,k,s) = f e st 5cj)(r,m,k,t)dt. 

Eventually, one obtains the following generalised eigenvalue problem, 

A 5<j)(r,m,k, s) = sL *5(j)(r,m,k,s), (47) 

where A and L* are 10M x 10M matrices, a few lines of which are subse¬ 
quently substituted for the implementation of the (no-slip) boundary condi¬ 
tions on the velocity. 

Let s* be the maximal growth rate, i.e., the real part of the maximal 
eigenvalue!! of Eg. 1471 over all possible wavenumbers m and k. Then the base 
flow is (linearly) stable if, and only if, s* < 0. 
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3 Due to the discretisation, some spurious eigenvalues may pop up in the generalised 
eigenvalue problem (Eq. 14711 , but they can easily be eliminated because, unlike their phys¬ 
ical counterparts, they vary with the number of discretisation points. 
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When the shear-thinning parameter a vanishes, these linear equations 
reduce to those derived by Avgousti and Beris |34| for an Oldro yd- B fluid; 
as a matter of fact, we have spotted slight differences with Ref. [34l |. which 
we believe are typos in that publication. 

These equations involve three non-dimensional quantities: the shear- 
thinning parameter a , the “bare Weissenberg” number t q 7 *, and the relative 
Newtonian viscosity ^ s /no, on top of the (fixed) relative curvature e. 

5-4■ Linear stability analysis: general trends and application to dense sus¬ 
pensions 

5.4.1. Consistency of the algorithm for a —> 0, ry s —> 0 and influence of the 
shear rate 

To check the validity of our algorithm, we have first set the shear-thinning 
parameter a and the Newtonian viscosity rj s to values close to zero. A con¬ 
ventional Upper Convected Maxwell (UCM) model should then be recovered, 
in which case, following Ref. [30|, some axisymnretric modes (at m = 0) be¬ 
come unstable at sufficiently large applied shear rates. This is indeed the 
case in our simulations, as soon as the shear rate grows larger than a crit¬ 
ical shear rate comparable to that reported in Ref. j3C)i |. However, we also 
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(a) s* vs. applied shear rate 7 for vanishing Newto¬ 
nian viscosity {^‘/vo = 10 -6 ) (UCM model). 


Vs/Vo 

(b) s* vs. r/ s , at fixed shear rate (7 = 20 ). 


Figure 5: Maximal growth rate s * (over all modes) in the limit of a non- 
shear-thinning fluid, a = 3 • 10 —4 , r Q = 1.4. At 7 = 20, the most unstable 
modes are at m ~ 1, k ~ 70. The dashed red line separates stable base flows 
(below the line) from unstable ones (above). 


observe that some non-axisymmetric modes, associated to small azimuthal 
wavenumbers m = 0(1), are even more unstable than the axisymmetric 
ones, in accordance with the literature |3J |. Figure [5] shows the increase of 
the growth rate s* of the most unstable mode with the shear rate 7 measured 
at the inner cylinder (or the “bare Weissenberg” number r a 7 ). 


5.4-2. Stabilising effect of the Newtonian viscosity 

Including a Newtonian contribution to the stress, via a finite value of the 
viscosity g s , tends to stabilise the flow, as shown in Fig. I5bl This stabilising 
role has already been reported in the literature on the Oldroyd-B model (i.e., 
for a = 0 ) 


32j; we find that it holds true for shear-thinning fluids, i.e., 
when a departs from zero (see Fig. [ 6 |). 


5-4-3. Stabilisation through shear-thinning 

For the values of a of interest here, shear-thinning strongly suppresses the 
instability. Figures 018] show that, for the model parameters corresponding 
to the colloidal suspension at 0 e g = 0.519, but with vanishing Newtonian 
viscosity, the base flow lies deep within the stable region, whereas, in the 
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Figure 6 : Dependence of the maximal growth rate s* on the Newtonian 
viscosity rj s , for a shear-thinning fluid with the following parameters: a = 
7.7, T a = 1.4, 7 = 20. (The growth rates vary only weakly with k and m for 
these parameters; only modes m ^ 250 are considered here). 

absence of shear-thinning, i.e., were a vanishingly small, a linear instability 
would develop, with a peak of instability around k ~ 70, m ~ 1. 

We should however mention the existence of a limited range of values of a , 
outside the experimentally relevant window for dense colloidal suspensions, 
namely, 0.3 < a < 3 and marked by red vertical bars on Fig. [71 which display 
unstable modes associated to abnormally large growth rates. These modes 
are located in a very different region of the ( k , m)-plane, namely k ss 0 and 
m 3 > 1. Although these perturbations have reasonable shapes in space, it 
is plausible that they are actually unphysical, but we do not know whether 
they are intrinsic in our WM-model or whether they arise because of artifacts 
in the (well established) numerical method. In the following, we concentrate 
on the experimentally relevant range of parameters. 

5 . 4 - 4 - Linear stability analysis of dense colloidal suspensions 

We now specifically consider the model parameters used to fit the rheol¬ 
ogy of dense colloidal suspensions, (see Tabled]). 

Consistently with the strongly stabilising effect of shear-thinning de¬ 
scribed in the previous section ( Section I5.4.3p . we have numerically ascer¬ 
tained the stability of the base flows corresponding to the suspension at 
cj) = 0.519 for rescaled shear rates 7 = 10~ 4 , 10~ 2 , 10 _1 , at <f> e ^ = 0.626 for 
7 = 5- 10~ 6 , 5 • 10 , 5 • 10~ 3 , 10 ~ 2 in the experimental range, as well as 

that with effective volume fraction (f> e g = 0.641 at 7 = 10” 8 , 10~ 6 , 10~ 4 . 
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Figure 7: Dependence of the maximal growth rate s* on the shear-thinning 
parameter a, for r is/rj 0 = 10~ 6 , t q = 1.4, 7 = 20. The arrow points to the 
typical values of a of interest here. See the text for the description of the 
red vertical bars. 



Figure 8 : Colour maps of the maximal growth rates associated with each pair 
of wavenumbers (m, k) for two distinct, but small, values of a (see captions), 
at 7 = 20, with negligible Newtonian viscosity ( r is/r] 0 = 10~ 6 ). Stable modes 
appear in dark blue. 
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The relative contribution of the Newtonian viscosity to the stabilisation of 
the flow echoes the contribution of the Newtonian stress to the total stress: 
at high densities and low shear rates, the effect of the Newtonian viscosity 
is negligible. Admittedly, the more acute shear-thinning displayed by the 
WM models fitted to the densest suspensions results in an enhanced stabil¬ 
ity compared to the experimental systems; however the fact that they lie 
so deep in the stable region suggests that the materials they model, albeit 
somewhat less shear-thinning, are also stable with respect to the considered 
perturbations. This is further supported by the observed stability of the 
system at cf> = 0.519, which is both less shear-thinning than the others and 
very well described by the WM model. 


5.5. Rationalisation with the Pakdel-McKinley criterion 

Before we conclude, it is interesting to note that the stabilising effects of 
shear-thinning and of the Newtonian contribution are in line with the visco¬ 
elastic stability criterion proposed by Pakdel and McKinley in 1996 BS 
(less general versions of the criterion can be found in earlier publications^|30|). 
On the basis of a dimensional analysis of generic visco-elastic constitutive 
equations, these authors introduced a dimensionless number, written V here, 
for inertialess visco-elastic instabilities in curved geometries in analogy to the 
classical Taylor number for inertial instabilities, and propounded the idea 
that the base flow is susceptible to an instability if this dimensional number 
exceeds a given threshold of order unity. The Pakdel-McKinley number reads 


V= l JL^± 

n x ’ 

where l p = vqt is the typical distance travelled by a material substructure 
( e.g ., a polymer chain) along the base-flow streamline while relaxing, 1Z is the 
radius of curvature of the streamline, N\ is the first normal-stress-difference 
and X is the shear stress. For a UCM model, the ratio ^ is simply the 
Weissenberg number and measures the “anisotropy of the normal forces”, 
as phrased by Morozov and van Saarloos [381 ]. But even with the present 


WM model, which is somewhat more complicated, the criterion based on V 
seems to capture the observed trends: a Newtonian contribution to the stress 
increases X without altering N\. thereby reducing V and stabilising the flow. 
On the other hand, an enhanced propensity to shear-thinning results in a 
decrease of N\ as (1 + crr Q 7 ) -2 , while X only decreases as (1 + ar a 7) -1 (for 
r] s — > 0 ); as a result, V is reduced, which also explains the enhanced stability 
of the flow. 
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The Newtonian stress contribution and the propensity to shear-thinning 
thus act in synergy, but, for the rheology of the very dense colloidal suspen¬ 
sions of interest here, our results indicate that the impact of shear-thinning 
quantitatively prevails over that of the former. 


6. Discussion and Outlook 

In conclusion, we have made use of, and partially extended, a general 
theoretical framework based on projections onto density modes, to describe 
the rheology of dense colloidal suspensions in the vicinity of the glass tran¬ 
sition. Contrary to the previous theoretical works, which focused on strictly 
homogeneous flows, special attention was paid to the fate of spatial inho¬ 
mogeneities. However, the intricacy of the formalism forced us to resort 
to particularly strong approximations before we could study the stability 
of the flow in curved (Taylor-Couette) geometry. At the expense of these 
approximations, constitutive equations falling in the White-Metzner class 
were obtained; the resulting model was shown to capture the shear-thinning 
properties of the material and to be reasonably consistent with experimental 
measurements of the linear rheology and steady-state rheology of the suspen¬ 
sions, although not in a strictly quantitative way. Eventually, we analysed 
the stability of the visco-elastic flow and brought evidence that, in the exper¬ 
imental range of shear rates, shear-thinning (and to a much lesser extent the 
Newtonian stress contribution) strongly stabilise the flow. This may explain 
why visco-elastic instabilities have been observed in a variety of visco-elastic 
fluids, but not in the dense suspensions under consideration here: the flow 
strains and destroys the microstructure of the material, so much so that, if a 
material volume is displaced by a perturbation, the memory of the stress that 
it carries, through its microstructure, is suppressed by the flow too quickly 
to allow the possibility of a feedback mechanism. 
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Appendix A. Operator formulation of the switch to the auxiliary 
frame 

In Section 12.21 the switch to the auxiliary frame at a specific point was 
presented as a change of coordinates. However, it may also be conducted 
with the operator formalism. Indeed, consider an auxiliary system (denoted 
by tildes) with the following Smoluchowski operator 

N 

&(T,t) = ^2[di + Fi(T) + v(ri,t)} ■ d u 
i— 1 

where 

v (r, t) = v soW (r, t) - u solv (r 0 (t), t ). (A.l) 

As in Section 12.21 r a (t) is the position (given by Eq. [8]) of the “material 
point” advected by the solvent flow field, with position r Q at t a . 

In the auxiliary system, an observable g(r,T\t) evolves as 

dtg(r,T-,t) = (tf(r,t)g{r,T;t). 


Evaluating an observable in this auxiliary system is therefore equivalent to 
evaluating it in the auxiliary (denoted by primes) of Section [2.21 g(r,T]t) 
and g'(r, T'; t) have the same time derivative fP(r, t) = fl^ / (r / , t ) 


r'=r 


and they coincide at t a , so they are equal at all times. 


r'=r 
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The formal solution of Eq. IA.1I is 


f 0 &( r ,s)ds~( r p) 


9 (r,T;t) = e_ 

f* fit (r ,s)ds 


g(r, r). 

Using the shorthand (t) for u solv (r a (t),t) ■ YZiLi di, we observe that A^(t) 
commutes with all fit(s). Indeed, 


= v 


N 

£ 

i =1 

solv 



N 

\di + -fi (T) + V Solv (ri, s) -di 

v soW (r D (t),t) Y djg 

3 = 1 


N N 


(r 0 (t),t) ■ Y Y \ d o d i9 + dj (Fi ■ dig) + djV solv (r i: s) dig 


i =1 j =1 


= A ] (t)n\s)g, 

where we have made use of djV solv (rj, s) = 0 and Y2iLi Fi = 0, hence 

•N N 

j r- 


YliLi djFi = 0- It follows that 


g(r,T-,t) = e= /o t ^ v (*-»W.-)-i:iLia ie Ji t n t (r,.)cfa^ (rjr)j 

= e Z f ° tdsvSOlv(ro(s) ’ s) ' i: ?= 1 di g( r , T; t) 
and, if g does not depend intrinsically on space, 


In particular, the fluctuation advection term emerges when the equation is 
differentiated with respect to time. 

dtg(r,T;t) = e Z^ dsv (r ° (s) ’ s) ' dr d t g(r, T; t) - v solv (r Q (t),t) -d r g(r,T-t). 


This equation agrees with Eq. [lTJ 
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